Introducción
El presente documento contiene un breve descripción del flujo de trabajo realizado para la clasificación de coberturas a partir de imágenes satelitales. La clasificación de coberturas es una ardua tarea y más aún cuando las clases son tan similares entre sí como en el caso de los cultivos que pertenecen a una misma temporada, por ejemplo, soja temprana y tardía. Esto, junto a la relativa escasa cantidad de datos disponibles convierten esta tarea en una verdadero desafío. Un desafío no solo técnico sino también creativo.
Los distintos scripts usados en esta competencia se encuentran alojados en: https://github.com/alessiobocco/DesafioAgTech2020
Pasos previos
Antes de comenzar es necesario cargar todos los paquetes que serán utilizados en el presente documento. El análisis se realizó en R y Python por lo que se deben cargar los documentos para que estos lenguajes interactúen.
# Instalar el paquete pacman que permite instalar y/o cargar los paquetes necesarios
if (!require("pacman")) install.packages("pacman", repos = 'http://cran.us.r-project.org')
# Instalar o cargar los paquetes necesarios
pacman::p_load("dplyr", "here", "fs", "devtools", "glue", "readr", "sf", "progress", "ggfortify", "xts", "ggplot2", "lubridate", "leaflet", "reticulate", "sp", "raster",
"viridis", 'cowplot', 'knitr', 'kableExtra', 'purrr', 'reticulate')Caso de estudio
El primer paso en el análisis consiste en la descripción del área en estudio y su caracterización. Para ello, se mostrarán los datos de entrada provistos por los organizadores de la competencia y relevados por la Bolsa de Comercio de Rosario y la Bolsa de Cereales de Buenos Aires.
# Datos de train
train <- read.csv("/Users/alessiobocco/Documents/Data_Science/Data_Challenges/DesafioAgTech2020/dataset/data_train.csv") %>%
dplyr::mutate(train = 1)
# Datos de test
test <- read.csv("/Users/alessiobocco/Documents/Data_Science/Data_Challenges/DesafioAgTech2020/dataset/data_test.csv") %>%
dplyr::mutate(train = 0)
# Labels
labels <- read.csv("/Users/alessiobocco/Documents/Data_Science/Data_Challenges/DesafioAgTech2020/dataset/Etiquetas.csv")
# Dataset completo
dataset <- train %>%
rbind(test) %>%
dplyr::left_join(., labels, by = c('Cultivo')) %>%
dplyr::mutate(lon = Longitud, lat = Latitud) %>%
sf::st_as_sf(., coords = c( 'Longitud', 'Latitud'),
crs = '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0') %>%
dplyr::mutate(Campania = if_else(Campania == "18/19", 2018, 2019))
# Area de estudio
area_estudio <- sf::st_read("/Users/alessiobocco/Documents/Data_Science/Data_Challenges/DesafioAgTech2020/dataset/Gral_Lopez",
layer = 'Gral_Lopez') %>%
sf::st_transform(., '+proj=tmerc +lat_0=-90 +lon_0=-60 +k=1 +x_0=5500000 +y_0=0 +ellps=intl +twogs84=-148,136,90,0,0,0,0 +units=m +no_defs') %>%
sf::st_buffer(., dist = 10000) %>%
sf::st_transform(., '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')## Reading layer `Gral_Lopez' from data source `/Users/alessiobocco/Documents/Data_Science/Data_Challenges/DesafioAgTech2020/dataset/Gral_Lopez' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 14 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -62.87973 ymin: -34.38361 xmax: -61.12538 ymax: -33.40216
## geographic CRS: WGS 84
El dataset completo que contiene la combinación del set de train y test se muestra a continuación:
Id: Integer - Es el identificador del dataset particular para la salida a campo
Cultivo: Character - Es el tipo de cultivo que se encuentra en el punto. Estos puntos han sido trasladados manualmente para encontrarse dentro del lote. Esta variable es el objetivo de la clasificación.
Longitud: Numeric - Coordenada angular.
Latitud: Numeric - Coordenada angular.
Elevacion: Numeric - Elevación del terreno.
Dataset: Character - Proovedor de la verdad de campo. BC: Bolsa de Cereales de Buenos Aires; BCR: Bolsa de Comercio de Rosario.
Campania: Character - Campaña de la cual ha sido obtenida la verdad de campo. Existen dos campañas: 18/19 y 19/20.
GlobalId: Integer - Es el identificador global del dataset.
Losa puntos se encuentran distribuidos en el Departamento de General López, Santa Fe (Argentina). En el siguiente mapa se muestra la distribución espacial de los puntos. Los puntos rojos corresponden a los datos de test, es decir, no tienen un cultivo asociado mientras que los azules corresponden a los datos de train y por lo tanto tienen toda la información sobre el mismo. Al hacer click en cada uno de los puntos se muestran las principales características de cada uno de ellos.
# Rampa de colores
pal <- colorFactor(palette = 'RdYlBu', domain = dataset$train)
leaflet::leaflet(data = dataset) %>%
leaflet::addTiles() %>%
addMapPane("polygons", zIndex = 410) %>%
addMapPane("points", zIndex = 420) %>%
leaflet::setView(lat = -34, lng = -62, zoom = 8) %>%
addProviderTiles("Esri.WorldImagery") %>%
leaflet::addCircleMarkers(lat = ~lat, lng = ~lon, radius = 3,
stroke = FALSE, fillOpacity = 1, opacity = 1,
color = ~pal(train),
popup = ~sprintf("<b>%s (%s)</b><br>Lat.: %.3f<br>Lon.: %.3f<br>Elev: %.0f m <br>Campania: %.0f",
Cultivo, Tipo, lat, lon, Elevacion, Campania),
options = leafletOptions(pane = "points")) %>%
leaflet::addPolygons(data = area_estudio, weight = 2, opacity = 1, color = "white", dashArray = "3",
fillOpacity = 0.1,
options = leafletOptions(pane = "polygons")) %>%
leaflet::addLegend(pal = pal, values = ~train, opacity = 1, title = NULL,
position = "bottomright")## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'
Las coordenadas de estas observaciones son el punto de partida para la extracción de la información satelital pero antes de eso es necesario la descarga y procesamiento de las imágenes para esta región.
Descarga de imágenes
Existen numerosas fuentes de información satelital pública, cada una con características propias, con fortalezas y debilidades. Para esta tarea se prefirió el uso de Sentinel 2 ya que combina una muy buena resolución espectral con un tiempo de revisita menor a otras alternativas, lo que ayuda con problemas como la cobertura nubosa porque hay mayores probabilidades de obtener una imagen sin nubes. Para la creación de los mosaicos mensuales se utilizó el paquete rgee que permite conectar el ambiente de R con Google Earth Engine y realizar gran parte del procesamiento en los servidores de Google. Con esta herramienta se generó una serie mensual de imágenes, en el caso que hubiese más de una imagen con un 10% de nubosidad los valores de cada píxel se agregaron utilizando la mediana. Cabe mencionar que para facilitar el procesamiento posterior y disminuir la carga de memoria, las imágenes se agregaron a una resolución de 30 metros. No se descargaron todas las bandas disponibles pero si las siguientes:
- B2: Azul (0.490 um)
- B3: Verde (0.560 um)
- B4: Rojo (0.665 um)
- B5: Vegetación (0.705)
- B6: Vegetación (0.740)
- B7: Vegetación (0.783)
- B8: NIR (0.842)
- B8A: Vegetación (0.865)
- B10: SWIR - Cirrus (1.375)
- B11: SWIR (1.610)
En la presente Figura se muestra un ejemplo de un índice espectral calculado a partir de la imagen Sentinel 2 correspondiente a noviembre de 2018.
Como se mencionó más arriba, el preprocesamiento en Google Earth Engine no incluyó la creación de máscaras de nubes. Por ello se realizaron máscaras empíricas basadas en las bandas 2, 10 y 11. Se determinaron umbrales para cada una de estas bandas para detectar la presencia de nubes y sus sombras y mediante operaciones de ráster se eliminaron los pixeles cubiertos.
Una vez descargadas las imágenes, una por cada mes del año de las Campañas 2018/2019 y 2019/2020, se extrajo la información correspondiente a cada uno de los puntos de los set de train y test mostrados anteriormente.
# Se cargan las bandas correspondientes a cada uno de los puntos
load("/Users/alessiobocco/Documents/Data_Science/Data_Challenges/DesafioAgTech2020/dataset/dataset_satelite_observado.RData")A continuación se visualizan los valores de las bandas extraídas de las imágenes de satélite
## [1] 351250 13
Debido a la presencia de nubes existen datos faltantes que deben ser imputados para no perder las observaciones que de por sí son escasas. La presencia de faltantes puede verse en la siguiente Figura.
La Figura anterior muestra la cantidad de datos válidos por mes para cada uno de los set, train y test. Una imagen con todos sus píxeles válidos para las bandas descargadas tiene un total de 8500 pixeles para el set de train y 5550 para el set de test. Valores menores indican que las máscaras de nubes han removido píxeles no validos por lo que deben ser imputados. Para el proceso de imputación se utilizó la librería missForest que utiliza árboles para la imputación de datos faltantes. Se trabajaron de manera independiente ambos set de datos para evitar la contaminación del set de test.
Al analizar la cantidad de datos disponibles se llegó a la conclusión de que podrían aumentarse la cantidad usada para entrenar el modelo al considerar los píxeles vecinos a cada punto observado. Para ello, se generó una grilla regular de 30 m de lado para toda el área en estudio. Luego sobre cada punto del set de train se crearon polígonos de Voronoi para asignarle una superficie a cada uno y también para evitar solapamientos entre puntos. Estos polígonos fueron luego recortados con un buffer de 120 m alrededor de cada punto observado para obtener así el área que efectivamente corresponde al mismo cultivo. Para evitar la contaminación de los datos complementarios de train por sobre los de test se eliminaron todos los puntos con una geometría identifica a los datos originales.
## [1] 23157 259
Al incluir los vecinos más próximos a los puntos de train es posible incrementar la cantidad de datos para entrenar el modelo en gran medida al pasar de alrededor de 800 datos a más de 20 mil. La extracción del set de datos complementario de las imágenes ráster debió ser paralelizado para aprovechar la presencia de más núcleos y así disminuir el tiempo de procesamiento, que de otra forma, tomaba varias horas. Al igual que para los puntos originales, la presencia de nubes genera datos faltantes. En este caso, por la cantidad de datos complementarios utilizar librerías como missForest implica un tiempo computacinal excesivo. Por este motivo es que los faltantes se rellenaron tomando la mediana de cada banda en función del Tipo de cobertura. Es decir, si el pixel faltante de noviembre de 2018 correspondía a soja, se tomó la mediana de todo los pixeles de esa imagen que tenían soja. Otra característica del dataset es la confección de las series temporales de índices. Primero se definió un año agrícola que es diferente al año calendario, es decir, el primer mes del año comienza en julio mientras que el último es en agosto. Los puntos muestreados solo tendrán datos en el año agrícola correspondiente a su campaña. Es decir, una soja sembrada en Octubre de 2019 sólo tendrá datos para el año agrícola 2019, mientras que para los demás años tendrán valor 0.
Indices espectrales
A partir de las bandas descargadas es posible calcular una variedad de índices espectrales que resumen las distintas bandas en único valor y tienen como objetivo describir determinadas propiedades de las coberturas como su grado de verdor, la humedad de suelo, el estado de la vegetación, etc. Los índices calculados a partir de las imágenes fueron los siguientes:
- NGRDI: Índice Normalizado de la Diferencia Verde-Roja: (B5 - B4)/(B5 + B4)
- TGI: Índice Triangular de Verdor: −0.5 × [190 × (B4 − B3) − 120 × (B4 − B2)]
- NDVI: Índice Normalizado de Diferencia de Vegetación: (B8 − B4) / (B8 + B4)
- GNDVI: Índice Normalizado de Diferencia de Verdor: (B8 – B3) / (B8 + B3)
- EVI: Índice de Vegetación Mejorado: 2.5 * (B8 – B4) / ((B8 + 6 * B4 – 7.5 * B2) + 1))
- AVI: Índice de Vegetación Avanzado: [B8 * (1 – B4)*(B8 – B4)]^1/3
- SAVI: Índice de Vegetación Ajustado por el Suelo: (B08 – B04) / (B08 + B04 + 0,428) * (1,428)
- NDWI: Índice Normalizado de Diferencia de Agua: (B3 – B8) / (B3 + B8)
- MSI: Índice de Estrés de Humedad: B11 / B8
- BSI: Índice de suelo Desnudo: (B11 + B4) – (B8 + B2) / (B11 + B4) + (B8 + B2)
- MCARI: Índice Modificado del Ratio de Absorción de Clorofila: ((B5 - B4) - 0.2 * (B5 - B3)) * (B5 / B4)
- RED: Índice del borde Rojo: (B7 / B5)^-1
- NDRE: Índice Norma;izado de la diferencia del borde Rojo: (B8A - B5)/(B8A + B5)
- NDII: Índice Normalizado de la Diferencia Infrarroja: (B8A - B11) / (B8A + B11)
Creación de features
Una vez que las series temporales de índices espectrales estaban listas y sin faltantes se crearon variables que resumían su comportamiento a lo largo del tiempo. Se calcularon variables para cada año agrícola en su totalidad y para cada una de las estaciones astronómicas (verano, otoño, invierno y primavera). Las medidas de resumen fueron la mediana, mínimo, máximo y desvío estándar. El resultado se observa a continuación.
## [1] 24058 61